Carga paquetes
library(sf)
library(spdep)
library(tmap)
library(INLA)Ejemplo 2: Modelación de casos de COVID
Motivación
Al comenzar la pandemia producida por COVID-19, surgió el interés de analizar la distribución espacial de los riesgos infección.
Ajustamos un modelo lineal generalizado con efectos aleatorios espaciales con INLA para estimar riesgo de aparición del virus en diferentes radios censales de la Ciudad de Córdoba, Argentina.
Datos
Se trabajará con una base de datos vectorial (.gpkg) que incluye polígonos de los radios censales de Córdoba y Gran Córdoba. Para cada radio censal se cuenta con información de
- Cantidad de habitantes en miles de personas (Poblacion)
- Nivel de fragmentación urbano (fragmentacion)
- Valor de la Tierra promedio
- Porcentaje de hogares con nececsidades básicas insatisfechas (HogaresNBI)
- Porcentaje de hogares con hacinamiento (HogariesHacinamiento)
- Porcentaje de hogares con Jefe/Jefa con al menos nivel de estudio universitario (HogariesJefe.Univ.)
- Porcentaje de hogares con computadoras (HogaresConComput.)
- Cantidad de habitantes por \(km^2\) (HabSup)
- Cantidad de bancos dentro del radio censal (Bancos)
- Cantidad de farmacias dentro del radio censal (Farmacias)
- Cantidad de centros de salud dentro del radio censal (EstSaludInt)
- Cantidad de bancos, farmacias y centros de salud del radio censal (Ban_Farm_Salud)
- Cantidad de casos detectados de COVID-19 (Casos)
Lectura base de datos
datos <- st_read("data/Base_07_07_radios.gpkg", quiet = TRUE)- Cálculo de número de casos esperados por radio censal
datos$E <- datos$Poblacion * sum(datos$Casos) / sum(datos$Poblacion)tm_shape(datos) +
tm_polygons(fill = 'E')
- Cálculo de cociente de infección estandarizada (Standardized Infection Ratio, SIR) por radio censal
datos$SIR <- datos$Casos / datos$Etm_shape(datos) +
tm_polygons(fill = 'SIR')
- Referencia para el efecto aleatorio y el término del error
datos$re_u <- 1:nrow(datos)
datos$re_v <- 1:nrow(datos)Ajuste del Modelo
Generación de lista con vecindarios
# Generación de lista con vecindarios
nb <- poly2nb(datos)
head(nb)
plot(st_geometry(datos), border = "grey", lwd = 0.5)
plot(
nb,
coordinates(as(datos, "Spatial")),
add = TRUE,
col = "blue",
points = FALSE,
lwd = 0.5
)
## [[1]]
## [1] 2 3 4 29 37 237 251 252 258 266 267 1498 1500
##
## [[2]]
## [1] 1 3 8 14 37
##
## [[3]]
## [1] 1 2 4 8
##
## [[4]]
## [1] 1 3 5 6 7 8 266 267
##
## [[5]]
## [1] 4 6 267
##
## [[6]]
## [1] 4 5 7 10 11 270# Generación de vecindarios para INLA
file_adj <- tempfile("map.adj1")
nb2INLA(file_adj, nb)
g <- inla.read.graph(filename = file_adj)# Ajuste del Modelo inflado en ceros
summary(datos)
formula <-
Casos ~
Fragmentacion +
HogaresNBI +
HogaresHacinamiento +
HogaresJefe.Univ. +
Bancos +
f(re_u, model = "besag", graph = g) +
f(re_v, model = "iid")
res_inflpoi <-
inla(
formula,
family = "zeroinflatedpoisson1",
data = datos,
E = E,
control.predictor = list(compute = TRUE)
)
summary(res_inflpoi)## Poblacion Fragmentacion ValorTierra HogaresNBI
## Min. : 7.0 Min. :1.000 Min. : 200 Min. : 0.0000
## 1st Qu.: 629.8 1st Qu.:2.000 1st Qu.: 2450 1st Qu.: 0.3755
## Median : 835.0 Median :4.000 Median : 6000 Median : 1.0855
## Mean : 873.2 Mean :3.367 Mean : 17074 Mean : 1.6939
## 3rd Qu.:1086.0 3rd Qu.:4.000 3rd Qu.: 12500 3rd Qu.: 2.4040
## Max. :3273.0 Max. :4.000 Max. :220000 Max. :20.6897
## HogaresHacinamiento HogaresJefe.Univ. HogaresConComput. HabSup
## Min. :0.0000 Min. : 0.0000 Min. : 0.00 Min. : 0.006
## 1st Qu.:0.1364 1st Qu.: 0.9802 1st Qu.:11.96 1st Qu.: 44.619
## Median :0.5051 Median : 3.3943 Median :18.05 Median : 73.769
## Mean :0.8617 Mean : 5.5653 Mean :20.02 Mean : 83.815
## 3rd Qu.:1.2639 3rd Qu.: 9.0126 3rd Qu.:25.26 3rd Qu.: 97.168
## Max. :6.6606 Max. :26.6055 Max. :57.17 Max. :606.024
## Bancos Farmacias EstSaludInt Ban_Farm_Salud
## Min. :0.00000 Min. :0.0000 Min. :0.0000 Min. :0.0000
## 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000
## Median :0.00000 Median :0.0000 Median :0.0000 Median :0.0000
## Mean :0.05301 Mean :0.2446 Mean :0.1536 Mean :0.4512
## 3rd Qu.:0.00000 3rd Qu.:0.0000 3rd Qu.:0.0000 3rd Qu.:1.0000
## Max. :8.00000 Max. :6.0000 Max. :6.0000 Max. :9.0000
## Casos geom E SIR
## Min. : 0.0000 MULTIPOLYGON :1660 Min. :0.002241 Min. : 0.000
## 1st Qu.: 0.0000 epsg:22174 : 0 1st Qu.:0.201595 1st Qu.: 0.000
## Median : 0.0000 +proj=tmer...: 0 Median :0.267299 Median : 0.000
## Mean : 0.2795 Mean :0.279518 Mean : 1.007
## 3rd Qu.: 0.0000 3rd Qu.:0.347649 3rd Qu.: 0.000
## Max. :36.0000 Max. :1.047749 Max. :89.967
## re_u re_v
## Min. : 1.0 Min. : 1.0
## 1st Qu.: 415.8 1st Qu.: 415.8
## Median : 830.5 Median : 830.5
## Mean : 830.5 Mean : 830.5
## 3rd Qu.:1245.2 3rd Qu.:1245.2
## Max. :1660.0 Max. :1660.0
##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts,
## ", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
## scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
## ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
## verbose, ", " lincomb = lincomb, selection = selection, control.compute
## = control.compute, ", " control.predictor = control.predictor,
## control.family = control.family, ", " control.inla = control.inla,
## control.fixed = control.fixed, ", " control.mode = control.mode,
## control.expert = control.expert, ", " control.hazard = control.hazard,
## control.lincomb = control.lincomb, ", " control.update =
## control.update, control.lp.scale = control.lp.scale, ", "
## control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
## ", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
## num.threads, ", " keep = keep, working.directory = working.directory,
## silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug =
## debug, .parent.frame = .parent.frame)" )
## Time used:
## Pre = 0.407, Running = 9.17, Post = 0.373, Total = 9.95
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -0.663 0.476 -1.603 -0.661 0.265 -0.661 0
## Fragmentacion 0.171 0.105 -0.035 0.171 0.377 0.171 0
## HogaresNBI -0.058 0.070 -0.196 -0.058 0.079 -0.058 0
## HogaresHacinamiento 0.120 0.137 -0.148 0.119 0.387 0.119 0
## HogaresJefe.Univ. 0.051 0.022 0.007 0.051 0.095 0.051 0
## Bancos 0.135 0.193 -0.244 0.135 0.513 0.135 0
##
## Random effects:
## Name Model
## re_u Besags ICAR model
## re_v IID model
##
## Model hyperparameters:
## mean sd
## zero-probability parameter for zero-inflated poisson_1 5.85e-01 5.00e-02
## Precision for re_u 3.20e-01 5.50e-02
## Precision for re_v 2.18e+04 1.32e+04
## 0.025quant 0.5quant
## zero-probability parameter for zero-inflated poisson_1 0.481 5.87e-01
## Precision for re_u 0.223 3.16e-01
## Precision for re_v 7188.279 1.84e+04
## 0.975quant mode
## zero-probability parameter for zero-inflated poisson_1 6.78e-01 5.94e-01
## Precision for re_u 4.41e-01 3.09e-01
## Precision for re_v 5.68e+04 1.35e+04
##
## Marginal log-Likelihood: -2310.35
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')Visualización Espacial
Visualización Interactiva
datos$RR <- res_inflpoi$summary.fitted.values[, "mean"]
datos$RR_LI <- res_inflpoi$summary.fitted.values[, "0.025quant"]
datos$RR_LS <- res_inflpoi$summary.fitted.values[, "0.975quant"]
popup_vars <- c(
"Fragmentación (Nivel)" = 'Fragmentacion',
"Valor Tierra Medio ($/m2)" = 'ValorTierra',
"Hogares NBI (%)" = 'HogaresNBI',
"Hogares Hacinamiento (%)" = 'HogaresHacinamiento',
"Hogares Jefe Univ.(%)" = 'HogaresJefe.Univ.',
"Bancos",
"Casos",
"Población" = 'Poblacion',
"E",
"SIR",
"RR_LI",
"RR",
"RR_LS"
)tmap_mode('view')
mapas <-
tm_basemap(c(Urbano = "OpenStreetMap", Satelite = "Esri.WorldImagery")) +
tm_shape(datos,
name = 'Casos') +
tm_polygons(
fill = 'Casos',
fill.scale = tm_scale_intervals(
style = "fixed",
breaks = c(0, 1, 2, 3, 5, 10, 15, 25, 36),
values = "brewer.yl_or_br"),
fill.legend = tm_legend(title = "Casos por Radio",
scientific = TRUE,
format = "f",
digits = 0
),
col = "gray50",
col_alpha = .5,
popup.vars = popup_vars
) +
tm_shape(datos,
name = 'Tasa de Infección - SIR') +
tm_polygons(
fill = "SIR",
fill.scale = tm_scale_intervals(
style = "fixed",
breaks = c(0, 5, 10, 15, 20, 25, 30, 50, 60, 90),
values = "brewer.yl_or_br"),
fill.legend = tm_legend(
title = 'SIR',
scientific = TRUE,
format = "f",
digits = 0
),
col = "gray50",
col_alpha = .5,
popup.vars = popup_vars
) +
tm_shape(datos,
name = 'Riesgo Relativo') +
tm_polygons(
fill = "RR",
fill.scale = tm_scale_intervals(
style = "fixed",
breaks = c(0.4, 1, 2, 4, 8, 10, 15, 20, 30, 50, 85, 96),
values = "brewer.yl_or_br"),
fill.legend = tm_legend(
title = "Riesgo Relativo",
scientific = TRUE,
format = "f",
digits = 1
),
col = "gray50",
col_alpha = .5,
popup.vars = popup_vars
)
mapasimg <-
"https://ig.conae.unc.edu.ar/wp-content/uploads/sites/68/2022/04/cropped-cromas-68-1.png"
map <-
tmap_leaflet(mapas) %>%
leafem::addLogo(img,
url = "https://ig.conae.unc.edu.ar/",
width = 136,
height = 50) %>%
leaflet::addMiniMap(position = "bottomleft",
width = 150,
height = 150)
map
# mapshot(map, paste0("Mapa_radios_", format(Sys.time(), "%d_%m_%Y"), ".html"))